%% 问题2：运动员起跳到落地动力学模型
clear; clc;

%% 输入参数（以附表1中运动员1为例）
m = 68;         % 质量 (kg)
g = 9.8;        % 重力加速度 (m/s²)
k = 4951.05;    % 蹦床垂直刚度系数 (N/m)
h0 = 2.0;       % 起跳最高点高度 (m)，可调整
c = 10;         % 空气阻力系数 (N·s/m)，伸展姿势时可增大（如15），蜷缩时减小（如8）

%% 阶段1：下落阶段求解（从h0到z=0）
% 定义位移函数 z(t)
z_fun = @(t) h0 + (m*g/c)*t - (m^2*g/c^2)*(1 - exp(-c*t/m));
% 寻找落地时间t1（z(t1)=0）
t1 = fzero(@(t) z_fun(t), 0, 2);  % 初始猜测区间[0,2]秒
v_impact = (m*g/c)*(1 - exp(-c*t1/m));  % 落地瞬间速度（向下为正）

%% 阶段2：接触阶段求解（蹦床压缩过程）
omega = sqrt(k/m);                % 角频率
t2 = linspace(0, 2*pi/omega, 500); % 模拟时间（2个周期）

% 解析解：压缩量delta(t)
delta = (m*g/k).*(1 - cos(omega.*t2)) + (v_impact/omega).*sin(omega.*t2);
% 弹力F(t) = k*delta(t)
F = k.*delta;

%% 求解最大压缩量和最大冲击力
[delta_max, idx] = max(delta);
F_max = F(idx);
fprintf('落地时间 t1 = %.4f 秒\n', t1);
fprintf('落地速度 v_impact = %.4f m/s\n', v_impact);
fprintf('最大压缩量 delta_max = %.4f m\n', delta_max);
fprintf('最大冲击力 F_max = %.2f N\n', F_max);

%% 可视化分析
figure('Position', [100 100 1000 600]);

% 子图1：下落阶段速度与位移
subplot(2,2,1);
t_fall = linspace(0, t1, 100);
v_fall = (m*g/c)*(1 - exp(-c*t_fall/m));
z_fall = z_fun(t_fall);
plot(t_fall, v_fall, 'r-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('速度 (m/s)');
title('下落阶段速度变化');
grid on;

subplot(2,2,2);
plot(t_fall, z_fall, 'b-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('高度 (m)');
title('下落阶段高度变化');
grid on;

% 子图2：接触阶段压缩量与冲击力
subplot(2,2,3);
plot(t2, delta, 'g-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('压缩量 (m)');
title('蹦床压缩量随时间变化');
grid on;

subplot(2,2,4);
plot(t2, F, 'm-', 'LineWidth', 1.5);
xlabel('时间 (s)'); ylabel('弹力 (N)');
title('蹦床冲击力随时间变化');
grid on;

%% 敏感性分析：不同起跳高度h0的影响（示例）
h0_list = [1.5, 2.0, 2.5]; % 不同起跳高度
F_max_list = zeros(1, length(h0_list));
for i = 1:length(h0_list)
    h0 = h0_list(i);
    z_fun = @(t) h0 + (m*g/c)*t - (m^2*g/c^2)*(1 - exp(-c*t/m));
    t1 = fzero(@(t) z_fun(t), 0, 2);
    v_impact = (m*g/c)*(1 - exp(-c*t1/m));
    delta_max = sqrt((m*g/k)^2 + (v_impact/omega)^2);
    F_max_list(i) = k*delta_max;
end

figure;
plot(h0_list, F_max_list, 'bo-', 'LineWidth', 1.5);
xlabel('起跳高度 h0 (m)'); ylabel('最大冲击力 F_max (N)');
title('起跳高度对冲击力的影响');
grid on;